#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
交通大数据分析脚本（论文级完整版）
功能：模拟交通数据 → 数据清洗 → 多维度分析 → 高质量可视化（含统计图表+学术标注）
输出：CSV数据文件 + PNG/PDF可视化图表 + 分析报告文本
"""

# -------------------------- 环境配置与依赖检查 --------------------------
import os
import sys
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime, timedelta
from scipy import stats

# 配置学术风格图表
plt.style.use('seaborn-whitegrid')  # 白色网格背景
sns.set_palette("colorblind")  # 颜色盲友好配色
plt.rcParams['font.sans-serif'] = ['SimHei']  # 支持中文显示
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 检查依赖库
REQUIRED_LIBS = {
    'pandas': '1.3.0',
    'numpy': '1.21.0',
    'matplotlib': '3.4.0',
    'seaborn': '0.11.0',
    'scipy': '1.7.0'
}

# -------------------------- 核心数据生成函数 --------------------------
def generate_traffic_dataset(output_dir="data", days=30, seed=42):
    """
    生成包含多维度的模拟交通数据集（含时间、天气、事故等影响因素）
    参数：
        output_dir: 数据存储目录（自动创建）
        days: 模拟天数（默认30天）
        seed: 随机种子（保证可复现）
    返回：
        data_path: 生成的CSV文件路径
    """
    # 创建输出目录
    os.makedirs(output_dir, exist_ok=True)
    data_path = os.path.join(output_dir, "traffic_data.csv")
    
    # 生成时间序列（30天，每小时一条记录）
    start_date = datetime(2024, 1, 1)
    timestamps = [start_date + timedelta(hours=i) for i in range(days*24)]
    np.random.seed(seed)  # 固定随机种子
    
    # 模拟天气类型（晴/雨/雪/雾）
    weather_types = ["晴", "雨", "雪", "雾"]
    weather_probs = [0.7, 0.15, 0.1, 0.05]  # 各天气概率
    
    # 模拟交通事故（0=无，1=有）
    accident_prob = 0.03  # 每小时事故发生概率
    
    # 构建数据列表
    data = []
    for ts in timestamps:
        hour = ts.hour
        day_of_week = ts.weekday()  # 0=周一, 6=周日
        is_weekend = 1 if day_of_week >= 5 else 0  # 周末标记
        
        # 模拟天气（带概率权重）
        weather = np.random.choice(weather_types, p=weather_probs)
        
        # 基础流量（工作日/周末不同模式）
        if is_weekend:
            base_flow = 400 + 30 * np.sin(hour * 0.2)  # 周末流量波动较小
        else:
            base_flow = 600 + 50 * np.sin(hour * 0.2)  # 工作日基础流量更高
        
        # 高峰时段叠加（工作日早/晚高峰）
        if not is_weekend and (7 <= hour <= 9 or 17 <= hour <= 19):
            peak_factor = np.random.normal(1.8, 0.2)  # 高峰时段流量放大1.8倍（均值±标准差）
            base_flow *= peak_factor
        
        # 天气影响（雨/雪/雾降低流量）
        weather_effect = {
            "晴": 1.0,    # 晴天无影响
            "雨": 0.85,   # 雨天流量降低15%
            "雪": 0.7,    # 雪天流量降低30%
            "雾": 0.6     # 雾天流量降低40%
        }[weather]
        base_flow *= weather_effect
        
        # 事故影响（事故时流量激增或骤降）
        has_accident = np.random.binomial(1, accident_prob)
        if has_accident:
            accident_effect = np.random.choice([-0.3, 0.5])  # 30%概率事故导致拥堵（-30%），50%概率事故后疏散（+50%）
            base_flow *= (1 + accident_effect)
        
        # 确保流量非负
        traffic_volume = max(0, round(base_flow))
        
        # 添加到数据列表
        data.append({
            "timestamp": ts,
            "hour": hour,
            "day_of_week": day_of_week,
            "is_weekend": is_weekend,
            "weather": weather,
            "has_accident": has_accident,
            "traffic_volume": traffic_volume
        })
    
    # 转换为DataFrame并保存
    df = pd.DataFrame(data)
    df.to_csv(data_path, index=False)
    print(f"✅ 模拟数据集已生成（含30天×24小时数据）: {os.path.abspath(data_path)}")
    return data_path

# -------------------------- 数据清洗与预处理 --------------------------
def clean_traffic_data(data_path):
    """
    清洗交通数据（处理缺失值、异常值、格式错误）
    参数：
        data_path: 原始数据路径
    返回：
        cleaned_df: 清洗后的DataFrame
    """
    try:
        df = pd.read_csv(data_path)
    except Exception as e:
        raise ValueError(f"读取数据失败: {str(e)}")
    
    # 检查关键列是否存在
    required_columns = ["timestamp", "hour", "day_of_week", "weather", "has_accident", "traffic_volume"]
    if not set(required_columns).issubset(df.columns):
        missing = set(required_columns) - set(df.columns)
        raise ValueError(f"数据缺少必要列: {missing}")
    
    # 转换时间戳格式（确保为datetime类型）
    df["timestamp"] = pd.to_datetime(df["timestamp"], errors="coerce")
    invalid_timestamps = df["timestamp"].isna()
    if invalid_timestamps.any():
        print(f"⚠️ 检测到{invalid_timestamps.sum()}条无效时间戳，已删除")
        df = df[~invalid_timestamps]
    
    # 处理流量异常值（基于IQR方法）
    q1 = df["traffic_volume"].quantile(0.25)
    q3 = df["traffic_volume"].quantile(0.75)
    iqr = q3 - q1
    lower_bound = q1 - 1.5 * iqr
    upper_bound = q3 + 1.5 * iqr
    outliers = (df["traffic_volume"] < lower_bound) | (df["traffic_volume"] > upper_bound)
    if outliers.any():
        print(f"⚠️ 检测到{outliers.sum()}条流量异常值（超出IQR范围），已删除")
        df = df[~outliers]
    
    # 填充天气缺失值（假设缺失时为"晴"）
    df["weather"] = df["weather"].fillna("晴")
    
    # 确保布尔列正确（has_accident应为0/1）
    df["has_accident"] = df["has_accident"].astype(int)
    
    print(f"✅ 数据清洗完成，原始数据{len(pd.read_csv(data_path))}条 → 清洗后{len(df)}条")
    return df

# -------------------------- 多维度分析函数 --------------------------
def analyze_hourly_pattern(df):
    """分析小时级流量规律"""
    hourly_stats = df.groupby("hour")["traffic_volume"].agg(
        mean_flow="mean",
        max_flow="max",
        min_flow="min",
        count="count"
    ).reset_index()
    
    # 计算高峰时段（平均流量前3小时）
    peak_hours = hourly_stats.sort_values("mean_flow", ascending=False).head(3)["hour"].tolist()
    return hourly_stats, peak_hours

def analyze_daily_pattern(df):
    """分析日级流量趋势（按工作日/周末）"""
    daily_df = df.groupby(["day_of_week", "is_weekend"])["traffic_volume"].agg(
        mean_flow="mean",
        median_flow="median",
        std_flow="std"
    ).reset_index()
    
    # 计算周内流量变化（周一到周日）
    weekday_order = [0, 1, 2, 3, 4, 5, 6]  # 周一到周日
    # 创建辅助列存储自定义顺序的索引
    daily_df['weekday_idx'] = daily_df['day_of_week'].map(lambda y: weekday_order.index(y))
    # 按辅助列排序（使用 by 参数）
    weekday_df = daily_df[daily_df["is_weekend"] == 0].sort_values(by='weekday_idx')
    # 删除辅助列（可选，若不需要保留）
    weekday_df = weekday_df.drop(columns=['weekday_idx'])
    weekend_df = daily_df[daily_df["is_weekend"] == 1]
   
    return daily_df, weekday_df, weekend_df

def analyze_weather_impact(df):
    """分析天气对流量的影响"""
    weather_stats = df.groupby("weather")["traffic_volume"].agg(
        mean_flow="mean",
        max_flow="max",
        sample_size="count"
    ).reset_index()
    
    # 计算天气与流量的相关性（Spearman秩相关）
    corr, p_value = stats.spearmanr(df["weather"].map({"晴":0, "雨":1, "雪":2, "雾":3}), df["traffic_volume"])
    
    return weather_stats, corr, p_value

def analyze_accident_impact(df):
    """分析事故对流量的影响"""
    accident_df = df.groupby("has_accident")["traffic_volume"].agg(
        mean_flow="mean",
        std_flow="std",
        sample_size="count"
    ).reset_index()
    
    # 事故前后流量变化（滑动窗口）
    df["prev_accident"] = df["has_accident"].shift(1)  # 前1小时是否事故
    df["next_accident"] = df["has_accident"].shift(-1)  # 后1小时是否事故
    
    return accident_df

# -------------------------- 可视化函数 --------------------------
def save_figure(fig, output_path, dpi=300, bbox_inches="tight"):
    """保存图表（支持PNG/PDF格式）"""
    os.makedirs(os.path.dirname(output_path), exist_ok=True)
    fig.savefig(output_path, dpi=dpi, bbox_inches=bbox_inches)
    print(f"✅ 图表已保存: {os.path.abspath(output_path)}")

def plot_hourly_pattern(hourly_stats, peak_hours, df, output_dir="results"):
    """绘制小时级流量规律图"""
    fig, ax = plt.subplots(figsize=(12, 6))
    
    # 绘制柱状图（平均流量）
    bars = ax.bar(hourly_stats["hour"], hourly_stats["mean_flow"], 
                  color="skyblue", alpha=0.7, label="平均流量")
    
    # 标记高峰时段
    peak_bars = [bar for bar, hour in zip(bars, hourly_stats["hour"]) if hour in peak_hours]
    for bar in peak_bars:
        bar.set_color("red")
        bar.set_alpha(0.9)
    
    # 添加数值标签
    for bar in bars:
        height = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2., height + 50,
                f"{height:.0f}", ha='center', va='bottom', fontsize=8)
    
    # 图表装饰
    ax.set_title("各小时平均交通流量（含高峰时段标注）", fontsize=14, pad=20)
    ax.set_xlabel("小时（0-23）", fontsize=12)
    ax.set_ylabel("平均车辆数", fontsize=12)
    ax.set_xticks(range(24))
    ax.grid(axis="y", alpha=0.3)
    ax.legend()
    
    # 保存图表
    output_path = os.path.join(output_dir, "hourly_flow_pattern.png")
    save_figure(fig, output_path)

def plot_daily_pattern(weekday_df, weekend_df, df, output_dir="results"):
    """绘制日级流量趋势图（工作日vs周末）"""
    fig, ax = plt.subplots(figsize=(12, 6))
    
    # 绘制工作日趋势
    ax.plot(weekday_df["day_of_week"], weekday_df["mean_flow"], 
            marker="o", linewidth=2, color="green", label="工作日")
    
    # 绘制周末趋势
    ax.plot(weekend_df["day_of_week"], weekend_df["mean_flow"], 
            marker="s", linewidth=2, color="orange", label="周末")
    
    # 图表装饰
    ax.set_title("周内每日平均交通流量（工作日vs周末）", fontsize=14, pad=20)
    ax.set_xlabel("星期（0=周一, 6=周日）", fontsize=12)
    ax.set_ylabel("平均车辆数", fontsize=12)
    ax.set_xticks(range(7))
    ax.set_xticklabels(["周一", "周二", "周三", "周四", "周五", "周六", "周日"])
    ax.grid(axis="y", alpha=0.3)
    ax.legend()
    
    # 保存图表
    output_path = os.path.join(output_dir, "daily_flow_pattern.png")
    save_figure(fig, output_path)

def plot_weather_impact(weather_stats, corr, p_value, df, output_dir="results"):
    """绘制天气对流量影响的箱线图"""
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    
    # 箱线图：不同天气的流量分布
    sns.boxplot(x="weather", y="traffic_volume", data=df, ax=ax1, palette="colorblind")
    ax1.set_title("不同天气条件下的交通流量分布", fontsize=14, pad=20)
    ax1.set_xlabel("天气类型", fontsize=12)
    ax1.set_ylabel("车辆数", fontsize=12)
    
    # 散点图：天气等级与平均流量的相关性
    weather_order = ["晴", "雨", "雪", "雾"]
    ax2.scatter(
        x=pd.Categorical(weather_stats["weather"], categories=weather_order, ordered=True).codes,
        y=weather_stats["mean_flow"],
        s=100,
        color="purple",
        edgecolor="black"
    )
    ax2.set_xticks(range(len(weather_order)))
    ax2.set_xticklabels(weather_order)
    ax2.set_title(f"天气等级与平均流量的相关性（ρ={corr:.2f}, p={p_value:.4f}）", fontsize=14, pad=20)
    ax2.set_xlabel("天气等级（0=晴, 3=雾）", fontsize=12)
    ax2.set_ylabel("平均车辆数", fontsize=12)
    
    # 保存图表
    output_path = os.path.join(output_dir, "weather_impact_analysis.png")
    save_figure(fig, output_path)

def plot_accident_impact(accident_df, df, output_dir="results"):
    """绘制事故对流量的影响图"""
    fig, ax = plt.subplots(figsize=(12, 6))
    
    # 绘制事故与无事故的流量对比
    ax.bar(accident_df["has_accident"], accident_df["mean_flow"], 
           yerr=accident_df["std_flow"], capsize=5, 
           color=["lightcoral", "lightgreen"], alpha=0.7)
    
    # 添加文本标注
    text_str = "".join([
        f"事故样本量: {accident_df.loc[1, 'sample_size']}",
        f"无事故样本量: {accident_df.loc[0, 'sample_size']}",
        f"事故时平均流量: {accident_df.loc[1, 'mean_flow']:.0f}",
        f"无事故时平均流量: {accident_df.loc[0, 'mean_flow']:.0f}"
    ])
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax.text(0.5, -0.15, text_str, transform=ax.transAxes, 
            fontsize=10, verticalalignment='top', bbox=props)
    
    # 图表装饰
    ax.set_title("事故对交通流量的影响（均值±标准差）", fontsize=14, pad=20)
    ax.set_xlabel("是否发生事故（0=无, 1=有）", fontsize=12)
    ax.set_ylabel("平均车辆数", fontsize=12)
    ax.set_xticks([0, 1])
    ax.grid(axis="y", alpha=0.3)
    
    # 保存图表
    output_path = os.path.join(output_dir, "accident_impact_analysis.png")
    save_figure(fig, output_path)

# -------------------------- 主流程函数 --------------------------
def main():
    """主函数：整合全流程"""
    try:
        # ---------------------- 步骤1：生成模拟数据 ----------------------
        print("===== 步骤1：生成模拟交通数据 =====")
        data_path = generate_traffic_dataset(output_dir="data", days=30)
        
        # ---------------------- 步骤2：数据清洗 ----------------------
        print("===== 步骤2：数据清洗与预处理 =====")
        cleaned_df = clean_traffic_data(data_path)
        
        # ---------------------- 步骤3：多维度分析 ----------------------
        print("===== 步骤3：多维度数据分析 =====")
        
        # 小时级分析
        hourly_stats, peak_hours = analyze_hourly_pattern(cleaned_df)
        print("▶ 小时级流量规律：")
        print(hourly_stats[["hour", "mean_flow", "max_flow"]].to_string(index=False))
        print(f"高峰时段（前3）: {peak_hours} 点")
        
        # 日级分析
        daily_df, weekday_df, weekend_df = analyze_daily_pattern(cleaned_df)
        print("▶ 日级流量趋势：")
        print(daily_df[["day_of_week", "is_weekend", "mean_flow"]].to_string(index=False))
        
        # 天气影响分析
        weather_stats, corr, p_value = analyze_weather_impact(cleaned_df)
        print("▶ 天气对流量的影响：")
        print(weather_stats.to_string(index=False))
        print(f"天气等级与流量的Spearman相关系数: ρ={corr:.2f} (p={p_value:.4f})")
        
        # 事故影响分析
        accident_df = analyze_accident_impact(cleaned_df)
        print("▶ 事故对流量的影响：")
        print(accident_df.to_string(index=False))
        
        # ---------------------- 步骤4：可视化 ----------------------
        print("===== 步骤4：生成可视化图表 =====")
        os.makedirs("results", exist_ok=True)  # 创建结果目录
        
        # 小时级流量图（传递cleaned_df）
        plot_hourly_pattern(hourly_stats, peak_hours, cleaned_df, output_dir="results")
        
        # 日级流量图（传递cleaned_df）
        plot_daily_pattern(weekday_df, weekend_df, cleaned_df, output_dir="results")
        
        # 天气影响图（传递cleaned_df）
        plot_weather_impact(weather_stats, corr, p_value, cleaned_df, output_dir="results")
        
        # 事故影响图（传递cleaned_df）
        plot_accident_impact(accident_df, cleaned_df, output_dir="results")
        
        # ---------------------- 步骤5：生成分析报告 ----------------------
        print("===== 步骤5：生成分析报告 =====")
        report_path = os.path.join("results", "analysis_report.txt")
        with open(report_path, "w", encoding="utf-8") as f:
            f.write("=== 交通数据分析报告 ===")
            f.write("1. 数据概况：")
            f.write(f"   - 原始数据量: {len(pd.read_csv(data_path))} 条")
            f.write(f"   - 清洗后数据量: {len(cleaned_df)} 条")
            f.write(f"   - 数据时间范围: {cleaned_df['timestamp'].min()} 至 {cleaned_df['timestamp'].max()}")
            
            f.write("2. 关键发现：")
            f.write(f"- 高峰时段集中在 {peak_hours} 点（工作日早/晚高峰特征明显）")
            f.write(f"- 周末平均流量比工作日低约 {int((weekday_df['mean_flow'].mean() - weekend_df['mean_flow'].mean())/weekday_df['mean_flow'].mean()*100)}%")
            f.write(f"- 天气对流量影响显著（ρ={corr:.2f}），雾天流量最低（比晴天低约40%）")
            f.write(f"- 事故发生时流量平均增加 {int((accident_df.loc[1, 'mean_flow'] - accident_df.loc[0, 'mean_flow'])/accident_df.loc[0, 'mean_flow']*100)}%（可能与事故后疏散有关）")
            f.write("3. 图表说明：")
            f.write("- results/hourly_flow_pattern.png：各小时平均流量（红色标注高峰时段）")
            f.write("- results/daily_flow_pattern.png：周内每日流量趋势（绿色=工作日，橙色=周末）")
            f.write("- results/weather_impact_analysis.png：不同天气的流量分布与相关性")
            f.write("- results/accident_impact_analysis.png：事故对流量的影响（均值±标准差）")
        
        print(f"✅ 分析报告已生成: {os.path.abspath(report_path)}")
        
        # ---------------------- 完成提示 ----------------------
        print("===== 全流程执行完成 =====")
        print("生成文件列表：")
        print("- 数据文件: data/traffic_data.csv")
        print("- 分析报告: results/analysis_report.txt")
        print("- 可视化图表: results/*.png（共4张）")
        print("建议使用这些文件撰写论文，图表可直接插入正文，数据可用于统计分析。")
    
    except Exception as e:
        print(f"===== 执行过程中发生错误: {str(e)} =====")
        print("请检查：")
        print("1. 输出目录是否有写入权限（尝试手动创建data/results目录）")
        print("2. Python环境是否安装所有依赖库（运行check_dependencies()函数检查）")
        print("3. 数据生成是否成功（查看data/traffic_data.csv是否存在）")

# -------------------------- 执行主函数 --------------------------
if __name__ == "__main__":
    main()